Establish environment

First step is to establish the environment for us to be able to process the FlowCam data, this is done by running the setup function. In this chunk we start by clearing my Global environment, this is optional and can be removed if not desired.

Most importantly the chunk establishes the folder which we are going to work from, as well as loading in the additional meta-data we require from a .csv file.

The minimum data included in the .csv file is sampling time from start of experiment, and dilution of sample before it was run in the FlowCam, as well as columns with the Flask-number, Mesodinium acclimation, prey-treatment and replicate indicator. Hence the file-name ‘Dil_time.csv’.

Load the package “flowPeaks”

if (!require("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
## Bioconductor version '3.17' is out-of-date; the current release version '3.18'
##   is available with R version '4.3'; see https://bioconductor.org/install
BiocManager::install("flowPeaks")
## Bioconductor version 3.17 (BiocManager 1.30.22), R 4.3.0 (2023-04-21)
## Warning: package(s) not installed when version(s) same as or greater than current; use
##   `force = TRUE` to re-install: 'flowPeaks'
require(flowPeaks)
## Loading required package: flowPeaks

Load special functions

The Meso.flowPeaks function is a wrapper for flowPeaks that tailors its use to FlowCam data with Mesodinium and cryptophytes. The function takes a subset of data with just one organism type (either Mesodinium or cryptophytes) and partitions it using the k-means algorithm in flowPeaks. It then uses breakpoints (that can be set by the user, or which are automated) to identify the plastid colours associated with each cluster identified by flowPeaks.

The function outputs an annotated dataset as well as summary statistics.

Meso.flowPeaks <- function(dataset,flask,day,create_plot = F,unknown.partition = 5.5, bluegreen.partition = 4.5){
  
  # Function last updated: 21 November 2023
  
  # Function inputs:
  #  dataset = csv file from the FlowCam
  #  flask = ID number for the flask
  #  day = experimental day
  #  create_plot argument (T = makes plots, F = automatic setting; doesn't make plots)
  #  unknown.partition = channel 1 partition for "unknown" cell types. Default is at 5.5
  #  bluegreen.partition = channel 2 partition for bluegreen plastids. Default is at 4.5
  
  
  # Step 1: Partition dataset to just fluorescent particles
  dataset.fluo <- dataset[dataset$Ch1.Peak>0,]
  dataset.fluo <- dataset.fluo[dataset.fluo$Ch2.Peak>0,]
  
  # Step 2: Use flowPeaks package to fit and visualize the number of clusters of points
  invisible(capture.output(fp <- suppressMessages(flowPeaks(dataset.fluo[,c("Ch1.Peak",'Ch2.Peak')],tol=.1,h0=.4))))
      # The "capture.output" function captures the printed output of the flowPeaks function
      # The "invisible" function suppresses that output so it doesn't automatically print.
      # The "suppressMessages" function suppresses the messages output of flowPeaks

  # Step 3: Extract cluster assignments
  dataset.fluo$peak <- fp$peaks.cluster

  # Step 4: Assign clusters to colors (or unknown)
  fp.data <- as.data.frame(summary(fp))
  fp.data$plastid.color <- 'red' # Presume plastids are red in color
  for(i in 1:dim(fp.data)[1]){
    if(fp.data$Ch1.Peak.center[i] < unknown.partition){ # If they fall below the unknown threshold, assign them "black" color for unknown
      fp.data$plastid.color[i] <- 'black'
    } else{ #If not unknown, then assign them "blue" color for bluegreen plastid type
    if(fp.data$Ch2.Peak.center[i] < bluegreen.partition){
      fp.data$plastid.color[i] <- 'blue'
    }
  }}
  dataset.fluo$peak.Ch1 <- fp.data$Ch1.Peak.center[match(dataset.fluo$peak,fp.data$cluster.id)]
  dataset.fluo$peak.Ch2 <- fp.data$Ch2.Peak.center[match(dataset.fluo$peak,fp.data$cluster.id)]
  dataset.fluo$plastidcolor <- fp.data$plastid.color[match(dataset.fluo$peak,fp.data$cluster.id)]
  
  # Step 5: Map fluorescent data back onto the dataset for export.
  dataset.analyzed <- dataset
  dataset.analyzed$plastidcolor <- dataset.fluo$plastidcolor[match(dataset.analyzed$Capture.ID,dataset.fluo$Capture.ID)]
  if(dim(dataset.analyzed[is.na(dataset.analyzed$plastidcolor),])[1]>0){
  dataset.analyzed[is.na(dataset.analyzed$plastidcolor),]$plastidcolor <- 'nonfluorescent'}
  dataset.analyzed$plastidcolor <- factor(dataset.analyzed$plastidcolor,levels=c('red','blue','black','nonfluorescent'))
  
  # Step 6: Calculate summary statistics
  counts <- as.data.frame(t(as.matrix(table(dataset.analyzed$plastidcolor))))
  counts$tot.red <- counts$red + (counts$black + counts$nonfluorescent)*counts$red/(counts$red+counts$blue)
  counts$tot.blue <- counts$blue + (counts$black + counts$nonfluorescent)*counts$blue/(counts$red+counts$blue)
  counts$tot.meso <- counts$tot.blue+counts$tot.red
  Ch1.peaks <- as.data.frame(t(as.matrix(tapply(dataset.analyzed$Ch1.Peak,dataset.analyzed$plastidcolor,FUN=mean))))
  Ch2.peaks <- as.data.frame(t(as.matrix(tapply(dataset.analyzed$Ch2.Peak,dataset.analyzed$plastidcolor,FUN=mean))))
  Ch1.peaks.sd <- as.data.frame(t(as.matrix(tapply(dataset.analyzed$Ch1.Peak,dataset.analyzed$plastidcolor,FUN=sd))))
  Ch2.peaks.sd <- as.data.frame(t(as.matrix(tapply(dataset.analyzed$Ch2.Peak,dataset.analyzed$plastidcolor,FUN=sd))))
  
  # Step 7: Make plots if desired
  if(create_plot==T){
    
    # 1. Baseline plot of all data
    par(mfrow=c(2,2))
    par(mar=c(2,4,2.5,0.5))
    plot(dataset.fluo$Ch1.Peak,dataset.fluo$Ch2.Peak,ylab='Channel 2 peak (phycoerythrin)')
    text(min(dataset.fluo$Ch1.Peak),max(dataset.fluo$Ch2.Peak),dataset.fluo$Class[1],pos=4)
    text(min(dataset.fluo$Ch1.Peak),max(dataset.fluo$Ch2.Peak)-.15,paste(flask,', Day ',day),pos=4)
    
    # 2. Plot of flowpeaks package partitions
    par(mar=c(2,2,2.5,2.5))
    plot(fp)
    
    # 3. Plot of data with cluster assignments
    par(mar=c(4,4,.5,.5))
    plot(dataset.fluo$Ch1.Peak,dataset.fluo$Ch2.Peak,col=c('red','blue','cyan','green')[as.factor(dataset.fluo$peak)],xlab='Channel 1 peak (chlorophyll)',ylab='Channel 2 peak (phycoerythrin)')
    
    # 4. Plot of data assigned by plastid type
    par(mar=c(4,2,.5,2.5))
    plot(dataset.fluo$Ch1.Peak,dataset.fluo$Ch2.Peak,col=dataset.fluo$plastidcolor,xlab='Channel 1 peak (chlorophyll)')
    
  }

  Res <- list('dataset.analyzed' = dataset.analyzed, 'counts' = counts, 'Ch1.peaks' = Ch1.peaks, 'Ch2.peaks' = Ch2.peaks,'Ch1.peaks.sd' = Ch1.peaks.sd, 'Ch2.peaks.sd' = Ch2.peaks.sd)
  Res
  
  
}

Analyze FlowCam raw data to generate spreadsheet summarizing results

This chunk loops over the specified folders to generate a spreadsheet summarizing the Mesodinium and cryptophyte abundances by plastid type and experimental conditions.

#Using the previously defined Folder_list, a for loop is created to access each folder one at a time, and within those folder each sample is then accessed one at a time.
#Each data set is being processed using the MRMC.partitioning function, which distinguished Red Chloroplast MC from blue-green Chloroplast MC. By using the meta-data stored in the file-name a new data-sheet is created with the desired data where each row is from a single sample from a single day.

create_plot_meso = T #If TRUE, the partitioning will output partitioned Ch1 vs Ch2 plots -- MESODINIUM VERSION
create_plot_prey = T #If TRUE, the partitioning will output partitioned Ch1 vs Ch2 plots -- PREY VERSION

unk.part.Meso = 5.5 # unknown partition for Mesodinium
bg.part.Meso = 4.5 # blue-green partition for Mesodinium
unk.part.prey = 0 # unknown partition for cryptophytes
bg.part.prey = 4.5 # blue-green partition for cryptophytes

#The first for-loop is created to open the folder of the individual sampling day in the order as they were defined (usually chronologically)
for(dir_idx in 1:length(Folder_list)){
  
  #The working directory is set to the folder containing one days worth of data.
  current_folder <- paste(root_folder,"/",Folder_list[dir_idx],sep = "")
  
  #Recording the current experimental day
  today <- expt.day[dir_idx]
  
  #here the columns to retrieve timestamp and dilution is defined.
  dil_name <- paste("Dil.",today,sep = "")
  time_name <- paste("TS.",today,sep = "")
  
  #The function list.files(current_folder), will create a list of files that are located in the folder, and by using this list of files, a new for-loop is created accessing each file consecutively.
  for(i in 1:length(list.files(current_folder))){
    
    #extracting meta data from the file name
    #Characters 1 to 7 is the Flask index
    name_flask <- substring(list.files(current_folder)[i],
                            first = 1,
                            last = 7)
    
    #assigning flask index as a number
    num_flask <- as.integer(substring(list.files(current_folder)[i],
                                      first = 6,
                                      last = 7))
    
    #Characters 9 to 12 is the Mesodinium pre-experiment plastid type
    Meso_designation <- substring(list.files(current_folder)[i],
                                  first = 9,
                                  last = 12)
    
    #Characters 14 through 18 indicate the prey being fed during the experiment
    prey_designation <- substring(list.files(current_folder)[i],
                                  first = 14,
                                  last = 18)
    
    #Character 20 is the replicate indicator
    expt_rep <- substring(list.files(current_folder)[i],
                          first = 20,
                          last = 20)
    
    #The dilution factor is extracted from the meta-data data-sheet in the root-folder
    flask_dil <- expt.meta[expt.meta$Flask == num_flask,
                           names(expt.meta) == dil_name]
    
    #Relative time from beginning of experiment is extracted
    flask_time <- expt.meta[expt.meta$Flask == num_flask,
                           names(expt.meta) == time_name]
    
    medium <- expt.meta[expt.meta$Flask == num_flask,"Medium"]
    
    #The .csv file is loaded here, and extra columns are added with the additional meta-data required for the data-processing
    file.name <- paste(current_folder,"/",list.files(current_folder)[i],sep = "")
    
    vec <- read.csv(file.name) %>% 
      mutate(Flask = name_flask,
             replicate = expt_rep,
             Meso = Meso_designation,
             Prey = prey_designation,
             Day = as.numeric(substring(expt.day[dir_idx],2)))
    
    #Because the flow-cam keeps counting the seconds when resetting the pump, a constant is multiplied to the Elapsed time to get the correct volume aspirated. This constant however was calculated from using the 0.5 ml pump, so a new value is probably needed for the 1 ml pump.
    Asp.vol <- (max(vec$Elapsed.Time*0.967)/60)*0.150
    
    #R-vector containing the raw data-set with both Mesodinium and Crypotphyte data
    vec.raw <- vec
    
    Meso_filter <- vec.raw$Class == "Mesodinium"
    
    #Creating an R-vector that only contains the data classified as Mesodinium
    vec <- vec.raw[Meso_filter,]
    
    #Creating an R-vector that only contains the data classified as Cryptophyte (or rather it only contains data that has not been classified as Meosdinium)
    vec2 <- vec.raw[!Meso_filter,]
    
    #On the first sample of each experimental day, an empty spreadsheet is created for both Prey and Mesodinium. This will be populated as the for-loop iterates through each sample.
    if(i == 1){
      FCDP_meso <- vec
      FCDP_prey <- vec2
      
      #Creating data-frame for partitioned data results
      All.Results <- as.data.frame(1:length(list.files(current_folder)))
      
      colnames(All.Results) <- c("Flask")     # Columns 1-8: Metadata
      All.Results$Day <- NaN
      All.Results$Time <- NaN                 
      All.Results$Dilution <- NaN             
      All.Results$Aspirated.vol <- NaN        
      All.Results$Meso <- NaN                 
      All.Results$Prey <- NaN                 
      All.Results$Replicate <- NaN            
      All.Results$Medium <- NaN               
      All.Results$M.Red <- NaN               # Columns 10-16; from Res$counts
      All.Results$M.Blue <- NaN              
      All.Results$M.Undefined <- NaN         
      All.Results$M.Nonfluorescent <- NaN
      All.Results$M.TotRed <- NaN
      All.Results$M.TotBlue <- NaN
      All.Results$TotalMeso <- NaN          
      All.Results$M.Red.meanCh1peak <- NaN               # Columns 17-20; from Res$Ch1.peaks
      All.Results$M.Blue.meanCh1peak <- NaN              
      All.Results$M.Undefined.meanCh1peak <- NaN         
      All.Results$M.Nonfluorescent.meanCh1peak <- NaN
      All.Results$M.Red.sdCh1peak <- NaN               # Columns 21-24; from Res$Ch1.sds
      All.Results$M.Blue.sdCh1peak <- NaN              
      All.Results$M.Undefined.sdCh1peak <- NaN         
      All.Results$M.Nonfluorescent.sdCh1peak <- NaN
      All.Results$M.Red.meanCh2peak <- NaN               # Columns 25-28; from Res$Ch2.peaks
      All.Results$M.Blue.meanCh2peak <- NaN              
      All.Results$M.Undefined.meanCh2peak <- NaN         
      All.Results$M.Nonfluorescent.meanCh2peak <- NaN  
      All.Results$M.Red.sdCh2peak <- NaN               # Columns 29-32; from Res$Ch2.sds
      All.Results$M.Blue.sdCh2peak <- NaN              
      All.Results$M.Undefined.sdCh2peak <- NaN         
      All.Results$M.Nonfluorescent.sdCh2peak <- NaN  
      
      #Adding columns for cryptophytes
      All.Results$C.Red <- NaN               # Columns 9-15; from Res$counts
      All.Results$C.Blue <- NaN              
      All.Results$C.Undefined <- NaN         
      All.Results$C.Nonfluorescent <- NaN
      All.Results$C.TotRed <- NaN
      All.Results$C.TotBlue <- NaN
      All.Results$TotalCrypto <- NaN          
      All.Results$C.Red.meanCh1peak <- NaN               # Columns 16-19; from Res$Ch1.peaks
      All.Results$C.Blue.meanCh1peak <- NaN              
      All.Results$C.Undefined.meanCh1peak <- NaN         
      All.Results$C.Nonfluorescent.meanCh1peak <- NaN
      All.Results$C.Red.sdCh1peak <- NaN               # Columns 20-23; from Res$Ch1.sds
      All.Results$C.Blue.sdCh1peak <- NaN              
      All.Results$C.Undefined.sdCh1peak <- NaN         
      All.Results$C.Nonfluorescent.sdCh1peak <- NaN
      All.Results$C.Red.meanCh2peak <- NaN               # Columns 24-27; from Res$Ch2.peaks
      All.Results$C.Blue.meanCh2peak <- NaN              
      All.Results$C.Undefined.meanCh2peak <- NaN         
      All.Results$C.Nonfluorescent.meanCh2peak <- NaN  
      All.Results$C.Red.sdCh2peak <- NaN               # Columns 28-31; from Res$Ch2.sds
      All.Results$C.Blue.sdCh2peak <- NaN              
      All.Results$C.Undefined.sdCh2peak <- NaN         
      All.Results$C.Nonfluorescent.sdCh2peak <- NaN  
      
    } else{
      FCDP_prey <- rbind(FCDP_prey,vec2)
      FCDP_meso <- rbind(FCDP_meso,vec)
    }
    
    # Add relevant data for this sample
    All.Results$Flask[i] <- name_flask
    All.Results$Day[i] <- today
    All.Results$Time[i] <- flask_time
    All.Results$Dilution[i] <- flask_dil
    All.Results$Aspirated.vol[i] <- round(Asp.vol,3)
    All.Results$Meso[i] <- Meso_designation
    All.Results$Prey[i] <- prey_designation
    All.Results$Replicate[i] <- expt_rep
    All.Results$Medium[i] <- medium
    
    #Running the partitioning function for Mesodinium

    if(dim(vec)[1] > 0){ # If there are Mesodinium in the dataset
      Res <- Meso.flowPeaks(vec,name_flask,today,create_plot = create_plot_meso,unknown.partition = unk.part.Meso, bluegreen.partition = bg.part.Meso)
      #Assigning partitioned data for Mesodinium
      All.Results[i,10:16] <- Res$counts
      All.Results[i,17:20] <- Res$Ch1.peaks
      All.Results[i,21:24] <- Res$Ch1.peaks.sd
      All.Results[i,25:28] <- Res$Ch2.peaks
      All.Results[i,29:32] <- Res$Ch2.peaks.sd
    }
    
    if(dim(vec2)[1] > 0){ # If there are cryptophytes in the dataset
      #Running flowPeaks for Prey
      Res.crypto <- Meso.flowPeaks(vec2,name_flask,today,create_plot = create_plot_prey,unknown.partition = unk.part.prey, bluegreen.partition = bg.part.prey)
      #Assigning partitioned data for Prey
      All.Results[i,33:39] <- Res.crypto$counts
      All.Results[i,40:43] <- Res.crypto$Ch1.peaks
      All.Results[i,44:47] <- Res.crypto$Ch1.peaks.sd
      All.Results[i,48:51] <- Res.crypto$Ch2.peaks
      All.Results[i,52:55] <- Res.crypto$Ch2.peaks.sd
    }
    
    if(i == length(list.files(current_folder))){
        FCDP_prey <- FCDP_prey %>% mutate(Partition = "Prey")
        FCDP_meso <- FCDP_meso %>% mutate(Partition = "Mesodinium")
        FCDP_all <- rbind(FCDP_meso,FCDP_prey)
      
      
      Res_name <- paste(CSV_folder,"/Summary_",expt.day[dir_idx],".csv",sep = "")
      FCDP_name <- paste(CSV_folder,"/FCDP_all_",expt.day[dir_idx],".csv",sep = "")
      
      write.csv(All.Results,Res_name,row.names = F)
      write.csv(FCDP_all,FCDP_name,row.names = F)
    }
  }
  
  if(dir_idx == 1){
    All.Results.Concat <- All.Results
  } else{
    All.Results.Concat <- rbind(All.Results.Concat, All.Results)
  }
  
  if(dir_idx == length(Folder_list)){
    All_res_name <- paste("All_Summary",".csv",sep='')
    write.csv(All.Results.Concat,All_res_name,row.names = F)
  }
  
}

Examine experimental outcomes

ARC <- All.Results.Concat # Shortening the variable name


# Something strange happened with the first Aspirated Volume measurement that I don't understand. But based on the cryptophyte numbers, the volume should be more similar to the other aspiration volumes than the 0.199 that the datasheet recorded.
ARC$Aspirated.vol[1] <- mean(ARC$Aspirated.vol[2:29])

# Calculating concentrations of cells (in cells per mL), normalizing to the aspirated volume and adjusting for dilution
ARC$M.pmL <- ARC$TotalMeso/ARC$Aspirated.vol*ARC$Dilu
ARC$M.pmL.red <- ARC$M.TotRed/ARC$Aspirated.vol*ARC$Dilu
ARC$M.pmL.blue <- ARC$M.TotBlue/ARC$Aspirated.vol*ARC$Dilu

ARC$C.pmL <- ARC$TotalCrypto/ARC$Aspirated.vol*ARC$Dilu
ARC$C.pmL.red <- ARC$C.TotRed/ARC$Aspirated.vol*ARC$Dilu
ARC$C.pmL.blue <- ARC$C.TotBlue/ARC$Aspirated.vol*ARC$Dilu

# Calculating proportion of blue-green cryptophytes
ARC$C.propblue <- ARC$C.pmL.blue/ARC$C.pmL

# Rounding day to whole number so that we can group by timepoint in downstream calculations
ARC$Day.Numeric <- round(ARC$Time,digits=0)

Did M. rubrum eat prey? To answer this, we can look at prey dynamics over the course of the experiment.

# Colors from Hiroshige and Cassatt2, Metbrewer, https://github.com/BlakeRMills/MetBrewer
TAcol <- rgb(230/255,98/255,84/255)
TAcol.bg <- rgb(230/255,98/255,84/255,.5)
HPcol <- rgb(114/255,187/255,213/255)
HPcol.bg <- rgb(114/255,187/255,213/255,.5)
Bothcol <- rgb(144/255,113/255,158/255)
Bothcol.bg <- rgb(144/255,113/255,158/255,.5)
par(mar=c(4,4,1,1),mfrow=c(1,3))

# Panel 1: All cryptophytes
plot(jitter(ARC$Day.Numeric,.2),ARC$C.pmL/1000,las=1,xlab='Time (days)',ylab='Cryptophytes (1000s of cells/mL)',pch=21,bg=c(Bothcol.bg,'white')[as.factor(ARC$Meso)],col=Bothcol,cex=1.5)
legend(x='topleft',inset=c(0.03,0.03),legend=c('Control','+ M. rubrum'),pch=21,pt.cex=1.5,pt.bg=c(Bothcol,'white'))

Summ.Prey <- summarySE(data=ARC, 'C.pmL', groupvars=c('Day.Numeric','Meso'))
arrows(Summ.Prey$Day.Numeric,Summ.Prey$C.pmL/1000+Summ.Prey$se/1000,Summ.Prey$Day.Numeric,Summ.Prey$C.pmL/1000-Summ.Prey$se/1000,code=3,angle=90,length=.02)
points(Summ.Prey$Day.Numeric,Summ.Prey$C.pmL/1000,pch=21,cex=1.5,bg=c(Bothcol,'white')[as.factor(Summ.Prey$Meso)])

# Panel 2: Teleaulax
plot(jitter(ARC$Day.Numeric,.2),ARC$C.pmL.red/1000,las=1,xlab='Time (days)',ylab='T. amphioxeia (1000s of cells/mL)',pch=21,bg=c(TAcol.bg,'white')[as.factor(ARC$Meso)],col=TAcol,cex=1.5)
legend(x='topleft',inset=c(0.03,0.03),legend=c('Control','+ M. rubrum'),pch=21,pt.cex=1.5,pt.bg=c(TAcol,'white'))

Summ.Prey <- summarySE(data=ARC, 'C.pmL.red', groupvars=c('Day.Numeric','Meso'))
arrows(Summ.Prey$Day.Numeric,Summ.Prey$C.pmL.red/1000+Summ.Prey$se/1000,Summ.Prey$Day.Numeric,Summ.Prey$C.pmL.red/1000-Summ.Prey$se/1000,code=3,angle=90,length=.02)
points(Summ.Prey$Day.Numeric,Summ.Prey$C.pmL.red/1000,pch=21,cex=1.5,bg=c(TAcol,'white')[as.factor(Summ.Prey$Meso)])

# Panel 1: All cryptophytes
plot(jitter(ARC$Day.Numeric,.2),ARC$C.pmL.blue/1000,las=1,xlab='Time (days)',ylab='H. pacifica (1000s of cells/mL)',pch=21,bg=c(HPcol.bg,'white')[as.factor(ARC$Meso)],col=HPcol,cex=1.5)
legend(x='topleft',inset=c(0.03,0.03),legend=c('Control','+ M. rubrum'),pch=21,pt.cex=1.5,pt.bg=c(HPcol,'white'))

Summ.Prey <- summarySE(data=ARC, 'C.pmL.blue', groupvars=c('Day.Numeric','Meso'))
arrows(Summ.Prey$Day.Numeric,Summ.Prey$C.pmL.blue/1000+Summ.Prey$se/1000,Summ.Prey$Day.Numeric,Summ.Prey$C.pmL.blue/1000-Summ.Prey$se/1000,code=3,angle=90,length=.02)
points(Summ.Prey$Day.Numeric,Summ.Prey$C.pmL.blue/1000,pch=21,cex=1.5,bg=c(HPcol,'white')[as.factor(Summ.Prey$Meso)])

We can see that in the presence of M. rubrum, the prey are really reduced. This suggests they’re being eaten by the ciliate!

Did M. rubrum grow? To answer this, we can look at ciliate dynamics over the course of the experiment.

par(mar=c(4,4,1,1))

ARC.Meso <- ARC[ARC$Meso=='MRTA',]

plot(jitter(ARC.Meso$Day.Numeric,.4),ARC.Meso$M.pmL/1000,las=1,xlab='Time (days)',ylab='M. rubrum (1000s of cells/mL)',pch=21,bg=c('white'),col='gray50',cex=1.5)
#legend(x='topleft',inset=c(0.03,0.03),legend=c('Control','+ M. rubrum'),pch=21,pt.cex=1.5,pt.bg=c(Bothcol,'white'))

Summ.Prey <- summarySE(data=ARC.Meso, 'M.pmL', groupvars=c('Day.Numeric'))
arrows(Summ.Prey$Day.Numeric,Summ.Prey$M.pmL/1000+Summ.Prey$se/1000,Summ.Prey$Day.Numeric,Summ.Prey$M.pmL/1000-Summ.Prey$se/1000,code=3,angle=90,length=.02)
points(Summ.Prey$Day.Numeric,Summ.Prey$M.pmL/1000,pch=21,cex=1.5,bg='black')

OK, great, so we know that M. rubrum were growing (and we’ll need to account for this in any calculations of grazing rate). However, without a prey-free control, we can’t estimate how much the presence of prey enhanced M. rubrum’s growth rate.

Preferential Grazing

The real goal of this experiment, however, was to detect preferential grazing (if present).

Did M. rubrum preferentially eat one prey type? If so, we’d expect the relative abundances of the prey to change over time.

par(mar=c(4,4,1,1))

plot(jitter(ARC$Day.Numeric,.2),ARC$C.propblue,las=1,xlab='Time (days)',ylab='Proportion H. pacifica',pch=21,bg=c(Bothcol.bg,'white')[as.factor(ARC$Meso)],col=Bothcol,cex=1.5)
legend(x='topright',inset=c(0.03,0.03),legend=c('Control','+ M. rubrum'),pch=21,pt.cex=1.5,pt.bg=c(Bothcol,'white'))

Summ.Prey <- summarySE(data=ARC, 'C.propblue', groupvars=c('Day.Numeric','Meso'))
arrows(Summ.Prey$Day.Numeric,Summ.Prey$C.propblue+Summ.Prey$se,Summ.Prey$Day.Numeric,Summ.Prey$C.propblue-Summ.Prey$se,code=3,angle=90,length=.02)
points(Summ.Prey$Day.Numeric,Summ.Prey$C.propblue,pch=21,cex=1.5,bg=c(Bothcol,'white')[as.factor(Summ.Prey$Meso)])

A couple takeaways from this: First, the ratio stays the same between control and + M. rubrum treatments. This suggests that the cells are being eaten proportional to abundance, rather than with some preference for, e.g., T. amphioxeia. Second, the ratio is never 50:50, even at the start. This suggests it might be hard to get H. pacifica inoculated in the correct ratio (we might need to give it a “boost” in our inoculation calculations).

Although these data suggest that M. rubrum is killing HP and TA in proportion to their abundance, this doesn’t necessarily tell us that M. rubrum is retaining HP plastids when it has the choice between these plastids and TA ones. We can look at the Channel 2 fluorescence (a proxy for phycoerythrin content) to guess at this, but absent controls for M. rubrum this won’t be definitive.

par(mar=c(4,4,1,1),mfrow=c(1,2))

Chlcol <- rgb(127/255,160/255,116/255)
Chlcol.bg <- rgb(127/255,160/255,116/255,.5)
ARC.Meso <- ARC[ARC$Meso=='MRTA',]

plot(jitter(ARC.Meso$Day.Numeric,.4),ARC.Meso$M.Red.meanCh1peak,las=1,xlab='Time (days)',ylab='M. rubrum Ch 1 (Chlorophyll)',pch=21,bg=c('white'),col=Chlcol.bg,cex=1.5)
#legend(x='topright',inset=c(0.03,0.03),legend=c('Channel 1 (chl-a)','Channel 2 (PE)'),pch=21,pt.cex=1.5,pt.bg=c('gray20',TAcol))

Summ.Prey <- summarySE(data=ARC.Meso, 'M.Red.meanCh1peak', groupvars=c('Day.Numeric'))
arrows(Summ.Prey$Day.Numeric,Summ.Prey$M.Red.meanCh1peak+Summ.Prey$se,Summ.Prey$Day.Numeric,Summ.Prey$M.Red.meanCh1peak-Summ.Prey$se,code=3,angle=90,length=.02)
points(Summ.Prey$Day.Numeric,Summ.Prey$M.Red.meanCh1peak,pch=21,cex=1.5,bg=Chlcol)



plot(jitter(ARC.Meso$Day.Numeric,.4),ARC.Meso$M.Red.meanCh2peak,las=1,xlab='Time (days)',ylab='M. rubrum Ch 2 (Phycoerythrin)',pch=21,bg=c('white'),col=TAcol.bg,cex=1.5)
#legend(x='topright',inset=c(0.03,0.03),legend=c('Channel 1 (chl-a)','Channel 2 (PE)'),pch=21,pt.cex=1.5,pt.bg=c('gray20',TAcol))

Summ.Prey <- summarySE(data=ARC.Meso, 'M.Red.meanCh2peak', groupvars=c('Day.Numeric'))
arrows(Summ.Prey$Day.Numeric,Summ.Prey$M.Red.meanCh2peak+Summ.Prey$se,Summ.Prey$Day.Numeric,Summ.Prey$M.Red.meanCh2peak-Summ.Prey$se,code=3,angle=90,length=.02)
points(Summ.Prey$Day.Numeric,Summ.Prey$M.Red.meanCh2peak,pch=21,cex=1.5,bg=TAcol)